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ABSTRACT 

When transiting their host stars, hot Jupiters absorb about 10% of the light in the wings of 
p I 1 the stellar Lyman-alpha emission line. The absorption occurs at wavelengths Doppler- shifted 

rjQ \ from line center by ±100 km/s — larger than the thermal speeds with which partially neutral, 

• • ~10 4 K hydrogen escapes from hot Jupiter atmospheres. It has been proposed that the ab- 

■ sorption arises from ~10 6 K hydrogen from the host stellar wind, made momentarily neutral 

by charge exchange with planetary H I. The ±100 km/s velocities would then be attributed 
to the typical velocity dispersions of protons in the stellar wind — as inferred from space- 
craft measurements of the Solar wind. To test this proposal, we perform 2D hydrodynamic 
simulations of colliding hot Jupiter and stellar winds, augmented by a chemistry module to 
compute the amount of hot neutral hydrogen produced by charge exchange. We observe the 
contact discontinuity where the two winds meet to be Kelvin-Helmholtz unstable. The Kelvin- 
Helmholtz instability mixes the two winds; in the mixing layer, charge exchange reactions 
establish, within tens of seconds, a chemical equilibrium in which the neutral fraction of hot 
stellar hydrogen equals the neutral fraction of cold planetary hydrogen (about 20%). In our 
simulations, enough hot neutral hydrogen is generated to reproduce the transit observations, 
and the amount of absorption converges with both spatial resolution and time. Our calcu- 
lations support the idea that charge transfer between colliding winds correctly explains the 
Lyman-alpha transit observations — modulo the effects of magnetic fields, which we do not 
model but which may suppress mixing. Other neglected effects include, in order of decreasing 
importance, rotational forces related to orbital motion, gravity, and stellar radiation pressure; 
we discuss quantitatively the errors introduced by our approximations. How hot stellar hydro- 
gen cools when it collides with cold planetary hydrogen is also considered; a more careful 
treatment of how the mixing layer thermally equilibrates might explain the recent detection 
of B aimer Ha absorption in transiting hot Jupiters. 

Key words: stars: winds, outflows - planets and satellites: atmospheres - line: formation - 
methods: numerical - ultraviolet: planetary systems 



1 INTRODUCTION 

Gas-laden planets lose mass to space when their upper atmo- 
spheres are heated by stellar ultraviolet (UV) radiation. Ubiqui- 
tous in the Solar System, thermally-driven outflows modify the 
comp ositions of their under lying atmospheres over geologic time 
(e.g., Weissman et IDEL999). Thanks to the Hubble Space Tele- 
scope (HST), escaping winds are now observed from extrasolar 
hot Jupiters: Jovian-sized planets orbiting at distances < 0.05 AU 
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from their host stars and bathed in intense ionizing fields. Spec- 
troscopy with HST reveals absorption depths of ~2-10% in var- 
ious resonance transitions (H I, O I, C II, Si III and Mg II) 
when the planet transits the star, implying ga s outflows that ex- 
tend for at l east several planetary radii (e.g.. IVid al-Madiar et al. 
2003 1 VM03:lvidal-Madiar et al .l2004lBen-Jaffell2007l:lBen-Jaffel 
20081: IVidal-Madiar et al.ll2008l:lLecavelier Pes Etangs et alfeOia 



Fossati et all bOld : iLinskv et all l201of) ~ Recent observations of 
HD 189733b also indicate temporal variations in H I Lyman- 
a absorption, possibly correlate d with stellar X-ray activity 
(iLecavelier des Etangs et al.l 120121) . These data promise to con- 
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strain the compositions of hot Jupit er atmospheres and the degrees 
to wh ich they are vertically mixed (ILiang et alJl2003l : IMoses etaD 

The HST observations of hot Jupiter winds are accompanied 
by theoretical studies that model planetary out fl ows starting from 
first principles (e.g ., lYellel 120041 : lYelld 12006, iTian et al J 120051 : 
iGarcia Munozll2007l : Murray-Clay, Chiang, & Murray 2009, M09). 
These ID hydrodynamic models generally agree that hot Jupiters 
like HD 209458b and HD 189733b are emitting M ~ 10 10 - 
10 11 g/s in mostly hydrogen gas . Three -dim ensional models in - 
clude [Lecavelier des Etan gs et al ] d2004h and Ijaritz et al ] (l2005h . 
who emphasize the importance of tidal forces. 

Do the models agree with the observations? Linsky et 
al. (2010) find that their observations of C II absorption in HD 
209458b can be made consistent with modeled mass loss rates, as- 
suming the carbon abundance of the wind is not too different from 
solar. More comparisons between observation and theory would 
be welcome — particularly for hydrogen, the dominant component 
of the wind. But the observations of H I absorption have proven 
surprisingly difficult to interpret. On the one hand, the original 
measurements by VM03 indicate substantial (~10%) absorption at 
Doppler shifts of +100 km/s from the center of the H I Lyman-a 
line. On the other hand, theory (e.g., M09) indicates that plane- 
tary outflows, heated by photoionization to temperatures T < 10 4 
K, blow only at ~10 km/s. How can such slow planetary winds 
produce significant absorption at +100 km/s? Measurements of 
blueshifted velocities as large as -23 km/s in the case of HP 
18973 3b only accentuate this problem jLecavelier des Etangs et al.l 

iHolmstrom et al.l (120081 . H08) propose that the observed ener- 
getic neutral H atoms arise from charge exchange between plan- 
etary H I and protons from the incident stellar wind. In this in- 
terpretation, the +100 km/s velocities correspond to the thermal 
velocities of 10 6 K hydrogen from the star — hydrogen which is 
made neutral by electron-exchange with planetary H I. The situ- 
ation is analogous to that of the colliding winds of O star binaries 
dStevens etal.l ll992: Lamberts et ai]l201 li and references therein). 
The H I Lyman-Qf absorption arises from the contact discontinuity 
where the two winds meet, mix, and charge exchange to produce 
hot neutral hydrogen. 

The calculations of H08, and those of the follow-up study 
bv lEkenback etlD J201CL E10), are based on a Monte Carlo al- 
gorithm that tracks individual "meta-particles" of neutral hydrogen 
launched from the planet. The meta-particles collide and charge 
exchange with stellar wind protons outside a presumed planetary 
magnetosphere, which is modeled as an "obstacle" in the shape of 
a bow shock. Good agreement with the Ly-ar observations is ob- 
tained for a range of stellar and planetary wind parameters, and for 
a range of assumed obstacle sizes. 

In this work we further test the hypothesis of charge ex- 
change first explored by H08 and E10. Our methods are com- 
plementary: instead of adopting their kinetic approach, we solve 
the hydrodynamic equations. We do not prescribe any obstacle to 
deflect the stellar wind, but instead allow the planetary and stel- 
lar winds to meet and shape each other self-consistently via their 
ram and thermal pressures. Some aspects of our solution are not 
realistic — we ignore the Coriolis force, the centrifugal force, stel- 



lar tidal gravity, and magnetic fields Q Our goal is to develop a first- 
cut hydrodynamic-chemical model of the contact discontinuity be- 
tween the two winds where material mixes and charge exchanges. 
Simple and physically motivated scaling relations will be devel- 
oped between the amount of H I absorption and the properties of 
the stellar and planetary winds. 

The plan of this paper is as follows. In ^2] we describe our 
numerical methods, which involve augmenting our grid-based hy- 
drodynamics code to solve the chemical reactions of charge ex- 
change, and specifying special boundary conditions to launch the 
two winds. In |3]we present our results, including a direct compar- 
ison with the H I hy-a transit spectra of VM03, and a parameter 
study to elucidate how the absorption depth varies with stellar and 
planetary wind properties. A summary is given in §4\ together with 
an assessment of the shortcomings of our study and pointers toward 
future work. 



2 NUMERICAL METHODS 

In 32.11 we describe the hydrodynamics code used to simulate the 
colliding planetary and stellar winds. In 32.2l we detail the charge 
exchange reactions that were added to the code. In 32.31 we outline 
our post-processing procedure for computing the Lyman-Qf trans- 
mission spectrum. As a convenience to readers, in 32.4l we re-cap 
the differences between our treatment of colliding winds and that 
of E10/H08. 



2.1 Hydrodynamics: Code and Initial Conditions 

Our simulations are performed with HERACLES dGonzalez et al.l 
|2007|) F1 a grid-based code using a second-order Godunov scheme 
to solve the Euler equations: 

^+V-(pV) = 
ot 

^ + V-[pV®V + p/] = 

ot 

^ + V.[(E + p)V] = 
ot 

^+V-(p Xi \) = . fl) 

Ot 

Here p, V, /?, and E are the mass dens ity, velocity, pressure, an d 
total energy density, respectively (e.g., IClarke & Carsw ell 2003). 
The code tracks abundances of individual species: x t is the mass 
fraction of the z' th species of hydrogen, where i e {1, 2, 3, 4} to cover 
four possible combinations of ionization state (either neutral or ion- 
ized) and temperature (either "hot" because it is arises from the star 
or "cold" because it arises from the planet). The outer product is 
denoted 0, and / is the identity matrix. 

All our simulations are 2D Cartesian in the dimensions x (stel- 
locentric radius) and y (height above the planet's orbital plane). 
Equivalently the simulations may be regarded as 3D, but with no 
rotation and with a star and a planet that are infinite cylinders 



For rec ent explorations of s tar-planet interactions including magnetic 
forces, see Cohen et al 1 d2011allbh . These simulations do not resolve the 
mixing layer interface between the stellar and planetary winds. 
2 http://irfu.cea.fr/Projets/Site_heracles/index.html 
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oriented parallel to the z-axis. At fixed computational cost, two- 
dimensional simulations enjoy better spatial resolution than three- 
dimensional simulations and thus better resolve the fluid instabil- 
ities at the interface of the two winds. The standard box size is 
(L XJ L y ) = (40R P , 60R P ), where R p = 10 10 cm is the planet radius. 
The number of grid points ranges up to (N x ,N y ) = (6400,9600); 
see Table 1 . 

The star and its wind are modeled after the Sun and the So- 
lar wind. The stellar wind is injected through the left edge of the 
simulation box; the densities, velocities, and temperatures in the 
vertical column of cells at the box's left edge are fixed in time. 
Stellar wind properties as listed in Table[T]are given for a stellocen- 
tric distance r = ri aunch; * = 5R Q , near the box's left edge. Here the 
stellar wind density, temperature, sound speed, and flow speed are 
set to n* = 2.9 x 10 4 cm" 3 , T, = 10 6 K, c* = 129 km/s (computed 
for a mean molecular weight equal to half the proton mass, appro- 
priate for an - 100% i onized hydrogen plasm a), and y* = 130 
km/s Jsheelev et all l997; Ou emerais et al.lu007l : see also lLemairel 
1201 lh . respectively. Our stellar wind parameters are such that the 
implied spherically symmetric (3D) mass loss rate is 1 x 10 12 g/s or 
2 X 10" 14 M o /yr. 

Our stellar wind parameters are similar to those of the "slow" 
Solar wind in the Sun's equatorial plane. Compare our choices with 
those of Ekenback et al. (2011), who adopt a stellar wind speed of 
450 km/s. Their speed is closer to that of the "fast" Solar wind 
which emerges from coronal holes. At Solar minimum, the fast 
wind tends to be confined to large heliographic latitudes (polar re- 
gions), but at Solar maximum, the coronal holes migrate to lower 
lati tudes and the fas t wind can more readil y penetrate to the eclip- 
tic (iKohl et al.ll 19981 : iMcComas et alEool S. Bale 2012, personal 
communication). Evaporating hot Jupiters like HD 209458b and 
HD 189733b have orbit nor mals that are near l y aligned with th e 
spin axes of their host stars dWinn et al.ll2005l : Iwinn et al.ll2006h . 
Because such planets reside near their stellar equatorial planes, the 
slow equatorial Solar wind seems a better guide than the fast, more 
polar wind; nevertheless, as noted above, the fast wind is known to 
extend to low latitudes, and the speeds and densities of both winds 
vary by factors of order-unity or more with time. 

The stellar wind velocity at the left boundary is not plane- 
parallel but points radially away from the central star (located out- 
side the box). The density, velocity, and temperature in each cell 
at the boundary are computed by assuming that the central star 
emits a spherical isothermal wind whose velocity grows linearly 
with stellocentric distance r and whose density decreases as 1/r 3 . 
These scalings, which are modeled after empirical Solar wind mea- 
surements (e.g., Sheeley et al. 1997) and which maintain constant 
mass loss rate with r, are used only to define the left-edge bound- 
ary conditions and are not used in the simulation domain. Outflow 
boundary conditions are applied at the top, bottom, and right edges 
of the box. 

As a final comment about our choice of stellar wind parame- 
ters, we note that they are valid for the left-edge boundary at r = 
Haunch,* = 5R Q — not for the planet's orbital radius of r = \0R Q . The 
left-edge boundary must be far enough away from the planet that 
the stellar wind properties at the boundary are well- approximated 
by their "free-stream" values in the absence of any planetary ob- 
stacle. We will see in |3]that the stellar wind slows considerably 
between r = 5R Q and r = \0R Q as a consequence of the oncom- 
ing planetary wind. This region of deceleration is absent from the 
models of H08 and E 10. 



A circle of radius di au nch,p = 4R P , centered at position (l x , l y ) = 
(30/? p , 30R P ) (where the origin is located at the bottom left corner 
of the domain), defines the boundary where the assumed isotropic 
and radial planetary wind is launched. The properties of our sim- 
ulated planetary wind, which are similar to those of the stan- 
dard isurjersomcjno^s ofHD20945 8b by Garcfa-Munoz (2007) 
and iMurrav-Cla v et al are listed in Table Q] and are con- 

stant in time along the circular boundary. The density and veloc- 
ity of the planetary wind at this boundary are such that if the 
wind were spherically symmetric, the mass loss rate would be 
1.6 x 10 11 g/s. This value lies within the range estimated from 
observations by Linsky et al. (2010) and from energetic con- 
sidera tions (e.g jLecavelier des Etangsll2007l ; lEhrenreich & Desertl 
1201 lh . N ote that 1 - /+ = 20% of the planetary wind at launch 
is neutral dMurrav-Clav et al.l l2009h and available for charge ex- 
change. This neutral fraction represents a balance between pho- 
toionizations by extreme UV radiation and gas advection of neu - 
trals at a planetary altitude of 4-5 R p dMurrav-Clav et al.|[2009h . 
The planetary and stellar winds are barely supersonic at launch 
(Mach numbers M p - 1.2 and M* = 1.01). 

Gravity is neglected, as are all rotational forces. The pres- 
sure p is related to the internal energy density e = E - pV 2 /2 
via p = (7 - l)e, where 7 = 1.01. That is, gas is assumed to be- 
have nearly isothermally. This isothermal assumption should not be 
taken to mean that the temperature is the same across the simulation 
domain; the temperature of the stellar wind at injection is r* = 10 6 
K, while that of the planetary wind is T p = 7000 K0 Rather, the 
two winds, as long as they remain unmixed, tend to maintain their 
respective temperatures as they rarefy and compress. In reality, the 
stellar wind can keep, in and of itself, a near-isothermal profile on 
length scales of interest to us because thermal conduction times (es- 
timated, e.g., using the Spitzer conductivity) are short compared to 
dynamical times. Treating the planetary wind as an isothermal flow 
is less well justified, as cooling by adiabatic expansion can be a sig- 
nificant portion of the energy budget (Garcfa-Munoz 2007; M09). 
Nevertheless the error incurred by assuming the planetary wind is 
isothermal is small for our standard model because the planetary 
wind hardly travels beyond its launch radius of 4R P before it en- 
counters a shock; thus rarefaction factors are small. Furthermore, 
as noted above, shock compression factors are modest because the 
speed of the planetary wind is only marginally supersonic. Where 
the stellar and planetary winds meet and mix, the code ascribes an 
intermediate temperature 10 4 K < T < 10 6 K. This temperature, as 
computed by HERACLES, is used only for the hydrodynamic evolu- 
tion; it is not used for computing either the charge exchange reac- 
tions ( 32.21 ) or the transmission spectrum ( 32.31 ). 

Each simulation is performed in two steps. First only the plan- 
etary wind is launched from its boundary and allowed to fill the en- 
tire domain for 2 x 10 5 seconds. Second the stellar wind is injected 
through the left side of the box, by suitable assignment of ghost 
cells. This two-step procedure was found to minimize transients. 
The simulations typically run for 2 x 10 6 s, which corresponds to 
~60 box-crossing times for the stellar wind in the horizontal direc- 
tion. 



3 In the standard model of M09, the temperature starts at ~10 4 K at 
a planetocentric radi us of 1.1R P — consistent with the observations by 
Balle ster etail 120071) — and cools to -3000 K at 4R p . The temperatures 
calculated bv lGarcfa Munozl <2007h at 4R P are 6000-7000 K. 
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Table 1. Parameters of the winds at launch, and of the simulation box 



Stellar Wind 



Planetary Wind 



Haunch,* - 5/?© 

= 2.9E4/cm 3 
r* = 1E6 K 
v* = 130 km/s 
tf = l 

c* = 129 km/s 
M* = 1.01 



^launch,/? — 4/?^ 

w p = 3.9E6/cm 3 
T p = 7000 K 
v p = 12 km/s 
J? =0.8 
c p = 10 km/s 



M 



1.2 



radial (jc) direction 



vertical (y) direction 



LJRp = 40 

N x = 50, 100, 200,400, 
800, 1600, 3200, 6400 



2.2 Charge Exchange 



L y /Rp = 60 

N y = 75, 150, 300, 600, 
1200, 2400, 4800, 9600 



Charge exchange consists of the following forward and reverse re- 
actions: 



H h +H °c 



Hi 



(2) 



Hot (subscript h) ionized (superscript +) hydrogen emitted by the 
star can collide with cold (subscript c) neutral (superscript 0) hy- 
drogen emitted by the planet, neutralizing the former and ionizing 
the latter while preserving their kinetic energies. The reverse reac- 
tion occurs with an identical rate coefficient p (units of cm 3 /s; j3 is 
the cross section multiplied by the relative velocity). 

We have added reaction ^ to HERACLES by integrating the 
following equations in every grid cell (we refer to this portion of 
the calculation as the "chemistry step"): 

d(n H x° c ) _ d(n H x + h ) 

dt dt 

d(n H x + c ) d{n H x + h ) 

dt dt 

d(n H x° h ) _ d(n H x + h ) 



dt 

x + h + x° h + x + c + x° c = 1 . 



(3) 



Here n H is the total hydrogen number density (regardless of ioniza- 
tion state or temperature), and x^} is a number fraction (equiv- 
alently a mass fraction because the only element treated in the 
simulation is hydrogen). The rate coefficient J3 = 4 x 10~ 8 cm 3 /s 
is calculated by combining the energy-dependent cross section of 
iLindsav & Stebbingsl d2005h with a Maxwellian distribution for the 
relative velocity between hydrogen atoms at the two (constant) tem- 
peratures 7* and T p . The finite-difference forms of equations (0 are 

+(n+l) _ +(n) _ 7 / 0(n+l) +(n+l) _ +(n+l) 0(n+l)\ 
X h X h ~ ° \ X h X c x h X c ) 



JKn+D . 



JKn) . 



0(n+l) _ 0(n) , 



+(n+l) 
+(«+!) 



+ X, 



+ X-, 



h 

+(n) 
h 

+(n) 



i 



(4) 



where the superscript (n) refers to the n th time step, b = j3n H At, and 
A£ is the integration time step of HERACLES. Because the righthand 



side of the first of these equations is evaluated at step (n + 1) instead 
of step (n), our scheme is implicit. The first equation combines with 
the others to yield 



+(n) 



+(n+l) 



i +b 



(5) 



from which the remaining number fractions at time step (n + 1) 
are derived. Because our solution is implicit, the dimensionless 
timestep b can exceed unity (as it does for our runs at low spatial 
resolution), and the system will still relax to its correct equilibrium. 
This chemical equilibrium is discussed further in 33.2.31 

Note that in contrast to H08 and E10, our calculations account 
for the reverse reaction + — > + H^. Accounting for the 
reverse reaction helps us to avoid overestimating the amount of hot 
neutral hydrogen. Our calculations of rP h are still overestimated, 
however, because we neglect thermal equilibration, i.e., cooling of 
hot hydrogen by collisions with cold hydrogen. In 34.1l we estimate 
the error incurred to be on the order of unity. 

Our calculations of the neutral fraction in the mixing layer do 
not explicitly account for photoionizations by Lyman continuum 
photons, radiative recombinations, or advection of neutral hydro- 
gen from the planetary wind — bu t these effects are already in- 
cluded by [Murray] Clav et all d2009h whose planetary wind param- 
eters we use; see 32.11 

2.3 Lyman-of Absorption 

The transmission spectrum in the Lyman-a line is post-processed, 
i.e., calculated after HERACLES has finished running. Both hot and 
cold neutral hydrogen (n° h and n° c ) contribute to the Lyman-Qf optical 
depth. It is assumed that the hot and cold neutral hydrogen do not 
thermally equilibrate (see 34. II where we question this assumption). 
Thus in computing the opacity due to hot hydrogen, we adopt a 
kinetic temperature of = 10 6 K, and in computing the opacity 
due to cold hydrogen we take T p = 7000 K. In each grid cell, 
the wavelength at line center is Doppler shifted according to the 
horizontal component of the bulk velocity (the observer is to the 
far right of the simulation box). Voigt line profiles are used with a 
damping constant (Einstein A c oefficient) equal to T = 6.365 x 10 8 
s _1 (e.g. JVerhamme et al .1120061) . 

For each wavelength A, the line-of- sight optical depth T A (y) 
is evaluated along each horizontal row of cells pointing to the star 
(lying between the white dashed lines in FiguresQ]and[2]). The total 
absorption is then computed as 



A(A) = <l-exp(-TVi)> 



(6) 



where () denotes a 1 -dimensional spatial average over y. Of course 
the star actually presents a circular disc, but because the simula- 
tion is only 2D, our simple ID average seems fair. The absorption 
profile A (A) can be computed for every snapshot (timestep) of the 
simulation. 



2.4 Differences Between This Work and E10/H08 

The main difference between our methods and those of E10/H08 
is that we numerically solve the equations of hydrodynamics in a 
2D geometry, whereas E10/H08 simulate collisions of hydrogen 
"meta-particles" in a more kinetic, 3D treatment. Neither we nor 
they compute magnetic forces explicitly. 

E10 include forces arising from the orbit of the planet about 



© 0000 RAS, MNRAS 000, 000-000 



Colliding Planetary and Stellar Winds 5 



the star, including the Coriolis force, the centrifugal force, and stel- 
lar tidal gravity. We do not. Our focus is on resolving mixing and 
charge exchange in the interface between the two winds. To that 
end, we solve for both the forward and reverse reactions of charge 
exchange (equations EHS, whereas E10/H08 solve only for the for- 
ward reaction. Our equations permit a chemical equilibrium to be 
established in the mixing layer; see 33.2.31 Furthermore, the struc- 
ture and geometry of the interaction region between the two winds 
are direct outcomes of our simulations, whereas the shape of the 
interface layer is imposed as a fixed "obstacle" in the simulations 
ofElO. 

Other differences include our treatments of the planetary and 
stellar winds. We account for both the neutral and ionized compo- 
nents of the planetary wind; E10 assume the planetary outflow is 
purely neutral. We draw our parameters of the stellar wind from 
those of the slow equatorial Solar wind, which blows at ~130 
km/s at a stellocentric distance of r = 5R Q (Sheeley et al. 1997; 
lOuemerais et al.ll2007h . E10 take the stellar wind to blow at 450 
km/s, while H08 take the stellar wind to blow at 50 km/s. Neither 
work accounts for how the stellar wind decelerates due to its in- 
teraction with the planetary wind, whereas in our simulations the 
deceleration zones are well-resolved. 

We will review again our simulation methods, and assess the 
severity of our approximations, in 34.11 



3 RESULTS 

Results for Lyman-a absorption by the mixing layer, including nu- 
merical convergence tests and a direct comparison with observa- 
tions, are given in 33.11 A parameter study is described in 33.21 

3.1 Absorption vs. Spatial Resolution and Time 

In Figures [T] and [2 we present results at our lowest (50 x 75) and 
near-highest (3200 x 4800) spatial resolutions, respectively. The 
simulations agree on the basic properties of the flow. The planetary 
wind is launched from the red circle and encounters a bow shock, 
visible in the left panels as a curved boundary separating orange 
(unshocked planetary wind) from red (shocked planetary wind). 
The radius of curvature of the planetary wind shock is roughly 
~6R P . Outside, the red region of thickness ~5R P contains shocked 
planetary wind. 

The stellar wind encounters a weak shock — visible as a near- 
vertical line separating dark blue from lighter blue in the left-hand 
panels of Figures [T]and[2] — at a distance of ~5R P from the left edge 
of the box. The shocked stellar wind is diverted around the planet 
by the pressure at the stagnation point where the two winds collide 
head on. 

We observe that both winds accelerate somewhat before they 
encounter shocks. For our standard model, the Mach numbers are 
M* < 1.3 and M p < 1.5 (for the parameter study simulations 
of 33.21 M p can grow up to 2-3). Density enhancements are thus 
modest — less than a factor of 2. 

The contact discontinuity between the stellar and planetary 
winds separates light blue from dark red in the left panels. It is lami- 
nar at low resolution but breaks up into turbulent Kelvin-Helmholtz 
rolls at high resolution (cf. Stone & Proga 2009 whose spatial res- 
olution was probably too low to detect the Kelvin-Helmholtz insta- 
bility). The middle panels plot the density of hot neutral hydrogen 



produced by charge exchange in the mixing layer. The "head" of 
the mixing layer, located near the stagnation point, spans only one 
or two grid cells in the low resolution simulation. The high resolu- 
tion simulation resolves much better the head of the mixing layer. 
Zoomed-in snapshots of the head will be presented in 33.21 

In Figure [3] the star-averaged absorption A at an equivalent 
Doppler velocity of + 100 km/s (redshifted away from the observer) 
is plotted against time for a range of spatial resolutions. From t = 
to 2 X 10 5 s, the planetary wind fills the simulation domain; the 
absorption quickly settles down to a value of ~2%. At these early 
times, only cold (T p = 7000 K) neutral hydrogen from the planet is 
available to absorb in Ly-a, and it is clearly insufficient to explain 
the absorption observed with HST. 

Starting at t = 2 x 10 5 s, the stellar wind is injected into the 
box. The absorption attains a first peak when the planetary and stel- 
lar winds reach a rough momentum balance and a mixing layer 
containing hot (7* = 10 6 K) neutral hydrogen is established. The 
height of the first peak decreases with each factor of 2 improvement 
in grid resolution until a resolution of 3200x4800 is reached. En- 
couragingly, all of the absorption values calculated in the various 
simulations converge at late times. 

The 3200x4800 run is the best behaved, with the absorp- 
tion holding steady at A « 9% for 10 6 s. Compared to all other 
simulations at lower resolution, the 3200x4800 run is the only 
one in which Kelvin-Helmholtz rolls appear (more on the Kelvin- 
Helmholtz instability in 33.2.21 ). 

We further tested the convergence of the 3200x4800 run by 
performing an even higher resolution simulation with 6400x9600 
grid cells. Because of the expense of such a simulation, the initial 
conditions of the 6400x9600 run were taken from the 3200x4800 
run at t = 10 6 s, and integrated forward for only 3 x 10 5 s (approx- 
imately 9 box crossing times for the stellar wind in the horizontal 
direction). The absorption values versus time for the 6400x9600 
run are overlaid in Figure [3] and are practically indistinguishable 
from those of the 3200x4800 run. Having thus satisfied ourselves 
that the 3200x4800 run yields numerically convergent results, we 
will utilize this grid resolution (0.0\25R P per grid cell length) for 
further experiments to understand the dependence of the absorption 
on input parameters, as described in 33.21 

Figure |4] plots the absorption spectrum for our standard 
3200x4800 simulation at t = 2 x 10 6 s. The absorption A is eval- 
uated at wavelengths offset from the central rest-frame wavelength 
of the Lyman-Qf transition by 9 Doppler- shift velocities Av. Ab- 
sorption at -50 km/s is stronger than at +50 km/s, a consequence 
of neutral, charge-exchanged hydrogen from the star accelerating 
from the stagnation point toward the observer. At larger velocities 
|Av| > 100 km/s, the spectrum is more nearly reflection- symmetric 
about Av = 0, because the broadening is purely thermal at r* = 10 6 
K. 

Figure \5\ displays the same information as in Figure |H but in 
the full context of the Hubble Space Telescope observations. The 
agreement between the modeled and observed in-transit spectra is 
encouraging. 

3.2 Scaling Relations for Absorption in the Mixing Layer 

To understand how absorption in the mixing layer depends on in- 
put parameters, we performed 3 additional simulations varying the 
launch properties n p , n*, v p and T p . The altered parameters are 
listed in Table [2] For all 3 simulations, the box size was main- 
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Figure 1. Snapshots of total density and velocity (left panel), density of hot neutral hydrogen (njj, middle panel), and temperature (right panel) of the 50x75 
simulation, with parameters listed in Table [T] Snapshots are taken at t = 2 x 10 6 s. The temperature map shown in the right panel is computed by HERACLES 
and used only to compute the hydrodynamic evolution; it is not used to compute the charge exchange reactions or the Lyman-a spectrum (see 32.2W2.3t . The 
two dashed white lines represent sightlines to the stellar limbs. 
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Figure 2. Same as Figure[T]and for the same simulation parameters but at a grid resolution of 3200x4800. Kelvin-Helmholtz rolls at the contact discontinuity 
first appear at this resolution. An even higher resolution of 6400x9600 yields the same star- averaged absorption; see Figure [3] 



tained at (L x ,L y ) = (40R P ,60R P ) and the grid resolution was 
(N X9 N y ) = (3200,4800). 

Note that our parameter study is not exhaustive. For example, 
none of the simulations listed in Table[2]varies the Mach number at 
launch of either the planetary or stellar wind. Actually we have per- 
formed simulations varying the Mach number of the stellar wind. 
These behave as we would expect — in particular, increasing M* in- 
creases the amount of absorption because of the increased compres- 
sion in the stellar shock. Nevertheless we elect not to include these 
extra simulations in our parameter study below. Magnetic fields, ne- 
glected by our simulations but certainly present in the stellar wind 
if not also the planetary wind, would stiffen the gas and prevent the 
kind of compression that we see when we raise M*. 

In the following subsections, we explain our numerical results 



on the properties of the mixing layer with order-of-magnitude scal- 
ing relations. The mixing layer's location is analyzed in 33.2.11 its 
thickness in 33.2.21 the densities of its constituent species in 33.2.31 
and the column density and absorptivity of hot neutral hydrogen in 
$32A\ 



3.2.1 Location of the mixing layer 

Along the line joining the planet to the star, the mixing layer — 
equivalently, the contact discontinuity — is located approximately 
where the two winds reach pressure balance: 

pM + cl)=p p (v 2 p +c 2 p ). (7) 
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Figure 3. Lyman-a absorption A (equation [6j, evaluated at a Doppler-shift velocity of +100 km/s from line center, versus time and spatial resolution. The 
absorption converges in time for all simulations, but only for grid resolutions of 3200x4800 or greater does a unique value for the absorption emerge. The 
3200x4800 simulation is also the lowest resolution run to resolve Kelvin-Helmholtz billows; see Figure [2] 
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Table 2. Launch parameters for the 3 additional 40R p x 60R p simulations 
at 3200 x 4800 resolution. The stellar parameters v* and are kept at their 
nominal values from Table [T] Note that none of the Mach numbers change. 





Nominal 


n p T 


n* i 


vpjp T 


11: 


2.9E4/cm 3 


2.9E4/cm 3 


9.7E3/cm 3 


2.9E4/cm 3 




3.9E6/cm 3 


1.2E7/cm 3 


3.9E7/cm 3 


3.9E6/cm 3 


y p 


12 km/s 


12 km/s 


12 km/s 


V3xl2 km/s 




7000 K 


7000 K 


7000 K 


21000 K 


n 


% = 0.11 


3x^o 


3xft 


3x^o 



°-200 -150 -100 -50 50 100 150 200 
Doppler-shift velocity Av [km/s] 

Figure 4. Lyman-*? absorption A versus Doppler-shift velocity Av from line 
center, evaluated for our standard 3200x4800 simulation at t = 2 x 10 6 s. 
Absorption at -50 km/s is stronger than at +50 km/s because of the bulk mo- 
tion of charge-exchanged neutral hydrogen streaming from the star toward 
the observer. The line wings at larger Doppler shifts are primarily thermally 
broadened at = 10 6 K. The absorption A « 9% at Av = ±100 km/s, in 
accord with HST observations; see Figure [5] 



In equation (0, quantities are evaluated near the mixing layer, not at 
launch. Note further that in equation and in equations to follow, 
we ignore the distinction between shocked and unshocked gas, as 
wind Mach numbers are near unity. Idealizing each wind velocity 
as constant, we substitute p p = M p l(2nv p d p ) andp* = M*/(27rv*<i*) 
into equation 0, as appropriate for the 2D circular winds in our 
simulations. Here d p measures distance from the planet, and d* 
measures distance from the star. Then the distance from the planet 
to the mixing layer — i.e., the approximate radius of curvature of the 
mixing layer — is given by 

d p = d^R (8) 
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Figure 5. Observed out-of-transit (highest blue curve) and observed in-transit (green curve) Lyman-Qf spectra, reproduced from Figure 2 of Vidal-Madjar et 
al. (2003). In the line "core" from -42 to +32 km/s, where interstellar absorption is too strong to extract a planetary transit signal, the flux is set to zero. Our 
theoretical in-transit spectrum (red curve) is computed by multiplying the observed out-of-transit spectrum by 1 - A, where A is plotted in Figured The 
agreement between the theoretical and observed in-transit spectra is good, supporting the idea that charge exchange between the stellar and planetary winds 
correctly explains the observed absorption at Doppler-shift velocities around ±100 km/s. 



where 

^ = M p tf p +<*)/v p 

For our standard m odel, ( R - %) ^ 0.11 . Note that for 3D spherical 
winds, d p = d* ^|<R dStevens et al.lll992r) . but this relation is not rel- 
evant for our 2D Cartesian simulations — we will use (|8} instead. 

The parameters in Table |2] were chosen to increase % by a 
factor of 3 compared to its value in our fiducial model. By equa- 
tion (|8]), when % = 3%), the mixing layer should be displaced 3x 
farther away from the planet compared to its location in our stan- 
dard model, assuming the star is far enough away that d* is essen- 
tially fixed at the star-planet separation. Figure [6] displays zoomed- 
in snapshots of the mixing layers for all simulations in Table [2] 
Looking at the /^-positions of the mixing layers, and recalling that 
the planet sits at l x = 30R P , we find that the layer is displaced 
(30 - 8)/(30 - 21) « 2 Ax farther away in the three new simulations 
as compared to the standard model. We consider this close enough 
to our expected factor of 3, given our neglect of the rather thick 
layers of shocked gas surrounding the mixing layer. 

Figure [7] shows density profiles for hot neutral hydrogen in 
the mixing layer for the three simulations plus our standard model. 



Densities are averaged over l y and plotted against l x . The fact that 
the mixing layers in the simulations having % = 3%) align in posi- 
tion confirms that K is the dimensionless parameter controlling the 
location of the mixing layer. 

3.2.2 Thickness of the mixing layer 

Figure [7] also indicates that the thickness of the mixing layer, L mix , 
varies when we change input parameters. Empirically, we find that 
the variations are consistent with the relation 

L mix ~0.1^(^J (10) 

where, as before, the distinction between shocked and unshocked 
gas densities is ignored. 

We can rationalize dTOt as follows. The timescale for a mode 
of wavelength /1 K h to grow exponentially by the linear Kelvin- 
Helmholtz instability (KHI) is given by 

t ~ ^kh Pp + P* ^ ^KH jPp_ ^ j ^ 

KH V* - V p 27T(p p p*) l l 2 27TV* y p* 

(e.g.,[chandrasekhar 196ll). We assume that the thickness of the 
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Figure 6. Zoomed-in snapshots of hot neutral hydrogen in the four simulations used to study how the properties of the mixing layer depend on input parameters; 
see Table |3 Snapshots are taken near t = 2 x 10 6 for the standard model, and near t = 1 x 10 6 s for the others. The bottom white dashed line is the line of sight 
to the lower stellar limb. The upper two red dashed lines bracket the "sampling interval" over which the hot neutral hydrogen density is vertically averaged 
to produce the density profiles shown in Figure [7] (the uppermost red dashed line is also the sightline to the upper stellar limb). In cases (b), (c), and (d), the 
mixing layer is located farther from the planet (centered at l x = 30R p ) than is the case in (a). In cases (b) and (c), the horizontal thickness of the mixing layer 
is greater than in cases (a) and (d). And in case (c), the density of hot neutral hydrogen is lowest. See 33.2l for explanations. 
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Figure 7. Density profiles of charge-exchanged hot neutral hydrogen in the 
mixing layer, averaged vertically (between the dashed red lines in Figure |6) 
and plotted against horizontal position. The planet is located to the right at 
l x = 30R P . The mixing layers of the three non-standard simulations are all 
displaced farther from the planet than in the standard model, a consequence 
of increasing the ratio % of the momentum carried by the planetary wind 
to that of the stellar wind ( 83.2. It . The thicknesses of the mixing layers 
as shown by the red and green curves are larger than those shown by the 
blue and cyan curves, a consequence of changing the growth rate for the 
Kelvin-Helmholtz instability ( 33.2.21 >. The two dashed lines are predictions 
of equation (\6\ based on considerations of chemical equilibrium; the blue, 
green, and cyan curves correctly intersect the upper dashed line, while the 
red curve correctly intersects the lower dashed line ( 33.2.31 . 



mixing layer saturates when a certain mode first becomes non- 
linear. Near saturation, the velocity perpendicular to the back- 
ground shear flow becomes comparable to the shear flow velocity: 
v± ~ v* - v p ~ v*. Thus when the mode becomes nonlinear, the 
mixing layer has thickness L mix ~ v ± t Kii ~ Wkh/2tt) ^p p /p*. This 
result matches dTOt , if we assume the initial disturbance that de- 



velops into the mixing layer has a characteristic lengthscale that is 
fixed at /Ikh ~ R P - Our description of mode saturation can only ap- 
ply to locations not too far downstream from the stagnation point; 
far away, the flows are too strongly perturbed to be described by 
the linear growth timescale dTTb . 



3.2.3 Density of hot neutral hydrogen in the mixing layer 

The density of hot neutral hydrogen in the mixing layer is set by 
chemical equilibrium. Suppose that within the layer, the total den- 
sity ftH,mix is approximately the average of the planetary wind den- 
sity and the stellar wind density: 



n p + n* 



(12) 



2 2 

The densities in (TT2l are those of shocked gas, but as is the case for 
all of 33.21 we ignore for simplicity the difference in density be- 
tween pre-shock and post-shock gas (see 33.2.11 ). Because n p » ft*, 
charge exchange hardly alters the ionization state of the shocked — 
and still cold — planetary wind. That is, the values of x° c and x" c do 
not change as the dense planetary wind mixes with the dilute stel- 
lar wind. In particular, the ratio x®/x+ is fixed at its initial value of 

a -/;>//; = 1/4. 

The timescale for charge exchange is (ftH,mix/>) _1 ~ 10 s, much 
shorter than the hours required for stellar-occulting gas to travel 
from the stagnation point to regions off the projected stellar limb. 
Thus nearly all of the gas seen in transit is driven quickly into chem- 
ical equilibrium, which from equation (0 demands that: 



l - 



f + 



f + 



(13) 



(14) 



In other words, in the mixing layer, the ionization fraction of stel- 
lar wind material quickly slaves itself to the ionization fraction of 
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planetary wind material. Now all of the hot hydrogen (both neutral 
and ionized) in the mixing layer originates from the stellar wind; 
from <TT2t , we have 



( x h + x h) n a,mix ~ 2~ • 
Combining (TUT) with (TT5T) yields 

n° h = x h n UMx ~ \{±-fp)n* 



(15) 



(16) 



Equation dT6t is approximately confirmed by our numerical results 
in Figure 13 the horizontal dashed lines predicted by dT6t roughly 
match the densities from our numerical simulations. 



3.2.4 Column density and absorptivity of hot neutral hydrogen 

Combining (TTTTb with (IT6t gives the total column density of hot 
neutral hydrogen: 



0.05(1 -f;)(n p n*f /2 R p . (17) 



For our standard model, N% 

During planetary transit, the hot absorbing gas that covers the 
face of the star is located near the stagnation point. As such, the 
bulk line-of-sight velocity of transiting gas is much less than its 
thermal velocity, which is of order 100 km/s. Assuming that the 
gas is only thermally broadened, and that the gas is optically thin at 
wavelengths Doppler- shifted from line center by velocities Av, we 
construct an approximate, semi-empirical formula for the absorp- 
tion: 



A(Av) ~ <o- line _ ctr exp[-m^(Av) 2 /2^r,] 

~o.i(lz£)( "> ) m ( 

\ 0.2 /Uxl0 6 cm- 3 / \ 



(18) 



vl/2 



3 x 10 4 cm- 



10 6 K\ 1 /exp[-m^(Av) 2 /2^] 



0.5 



(19) 



where cr line _ ctr = 6 x 10 _15 (10 6 K/T*) 1/2 cm 2 is the line-center cross 
section for the Lyman-ar transition, m H is the mass of the hydro- 
gen atom, and k is Boltzmann's constant. Strictly speaking, the 
quantities in equation d~j~9T) should be evaluated in the vicinity of 
the contact discontinuity, but we have instead normalized equation 
dT9l ) to the wind properties at launch (evaluated at di au nch, P = 4R P 
and r launch; * = 5R Q ). We have verified in our simulations that the 
launch properties differ only by factors of order unity from the val- 
ues at the contact discontinuity, and so equation (\9\ may be used to 
predict the absorption by inserting only the launch properties. The 
exponential in equation (\9\ is evaluated for nominal parameters 
Av = 100 km/s and T+ = 10 6 K. 

As a further check, we show in Figure [8]the absorption values 
A to which the four simulations converge. They compare well with 
the values predicted by (\9l using only the launch properties. 

Had we kept the dependence of the mixing layer properties on 
the stellar wind Mach number M*, equation fj~6l) would be modified 
such that n° h oc M\n* — where n* is the pre-shock (launch) density — 
in accord with the usual Rankine-Hugoniot jump condition that 
states that the density increases by the square of the Mach number 
across a plane-parallel isothermal shock. And equations dTTb— dT9t 
would be modified such that A oc oc M*. Indeed our numeri- 
cal simulations (not shown) confirm this linear dependence of A on 
M* . We mention this result only in passing because it is not likely to 
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Figure 8. Lyman-a absorption A versus time, evaluated at a Doppler-shift 
velocity of +100 km/s from line center, for our standard model plus three 
additional models with different input parameters as indicated in the leg- 
end (see also Table UJ. The colored jagged lines are the results from our 
numerical simulations. The dashed black lines are the predictions from our 
physically motivated scaling relation (19) ; the simulations converge fairly 
well to the predicted values. 



remain true once we account for the real-life magnetization of the 
stellar wind. Magnetic fields stiffen gas and reduce the dependence 
of A on M*. 



4 SUMMARY AND DISCUSSION 

Using a 2D numerical hydrodynamics code, we simulated the col- 
lisional interaction between two winds, one emanating from a hot 
Jupiter and the other from its host star. The winds were assumed 
for simplicity to be unmagnetized. Properties of the stellar wind 
were drawn directly from observations of the equatorial slow So- 
lar wind (Sheeley et al. 1997; Quemerais et al. 2007; Lemaire 
2011), while those of the planetary wind were taken from hydro- 
dynamic models of outflows powered by photoionization heating 
(Garcia-Munoz 2007; Murray-Clay et al. 2009). For our standard 
parameters, the mass loss rate of the star is M* = 2 x 10 _14 M o /yr 
= 10 12 g/s and the mass loss rate of the planet is M p = 1.6 x 10 11 
g/s = 2.7 x 10 _3 Mj/Gyr. At the relevant distances, each wind is 
marginally supersonic — the stellar wind blows at ~ 130-170 km/s 
(sonic Mach number M* < 1.3) and the planetary wind blows at 
~ 12-15 km/s (Mach number M p < 1.5). Thus shock compression 
is modest, even without additional stiffening of the gas by magnetic 
fields. 

A strong shear flow exists at the contact discontinuity between 
the two winds. At sufficiently high spatial resolution, we observed 
the interfacial flow to be disrupted by the Kelvin-Helmholtz insta- 
bility. The Kelvin-Helmholtz rolls mix cold, partially neutral plan- 
etary gas with hot, completely ionized stellar gas. Charge exchange 
in the mixing layer produces observable amounts of hot (10 6 K) 
neutral hydrogen. Upon impacting the planetary wind, the hot stel- 
lar wind acquires, within tens of seconds, a neutral component 
whose fractional density equals the neutral fraction of the plane- 
tary wind (about 1 - f p = 20%). Seen transiting against the star, 
hot neutral hydrogen in the mixing layer absorbs ~10% of the light 
in the thermally broadened wings of the stellar Lyman-of emission 
line, at Doppler shifts of ~100 km/s from line center. Just such a 
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transit signal has been observed with the Hubble Space Telescope 
(Vidal-Madjar et al. 2003). The +100 km/s velocities reflect the 
characteristic velocity dispersions of protons in the stellar wind — 
as inferred from in- situ spacec raft observations of the Solar wind 
(e.g., Figure 3 of lMarschlEzOoih . 

Our work supports the proposal by Holmstrom et al. (2008) 
and Ekenback et al. (2010) that charge exchange between the 
stellar and planetary winds is responsible for the Ly-a absorp- 
tion observed by HST. This sa me conclusion is reached by 
iLecavelier des Etangs et alj d2Q12h in the specific case of HD 
189733b. Our ability to reproduce the observations corroborates 
the first-principles calculations of hot Jupiter mass loss on which 
we have relied (e.g., Yelle 2004; Garcia-Munoz 2007; Murray-Clay 
et al. 2009, M09). Time variations in Ly-a absorption are expected 
both from the variable stellar wind — the Solar wind is notoriously 
gusty — and from the variable planetary wind, whose mass loss rate 
tracks the time-variable ultraviolet and X-ray stellar luminosity. 



4.1 Neglected Effects and Directions for Future Research 

Although the general idea of photoionization-powered planetary 
outflows exchanging charge with their host stellar winds seems 
correct, details remain uncertain. We list below some unresolved 
issues, and review the effects that our simulations have neglected, 
in order of decreasing concern. 

(i) Thermal equilibration in the mixing layer Our calculations over- 
estimate the amount of hot neutral hydrogen produced by charge 
exchange because they neglect thermal equilibration. A hot neutral 
hydrogen atom cools by colliding with cold gas, both ionized and 
neutral, from the planetary wind. The concern is that hot neutral gas 
cools before it transits off the face of the star. Starting from where 
the mixing layer is well-developed (say the lower red dashed line in 
Figure [6]), hot neutral gas is advected off the projected stellar limb 
in a time 



fadv ~ 2R P /v* ~ 2 x 10 3 s . 
By comparison, the cooling time is of order 



(20) 



^cool 



1 



500 



2 x 10 6 cm" 3 



100 km/s 



(21) 

where n + c is the density of cold ionized hydrogen in the mixing 
layer, v re i is the relative speed between hot and cold hydrogen, and 
cr is the H-H + cross section for slowing down fast h ydrogen, here 
taken to be the "viscosity" cross section calculated bv lSchultz et alj 
(2008)0 Our estimate of t coo \ in (1211 neglects cooling by neutral- 
neutral collisions, but we estimate the correction to be small, as 
n° c is lower than n+ by a factor of 1/(1 - f£) ~ 5, and the cross 
section for H-H collisions is generally not greater than for H-H + 
collisions (A. Glassg old 2012 , personal communicat ion; see also 
ISwenson etalll 19851 : note that lEkenback et all (1201 Oh take the rel- 
evant H-H cross section to be 10~ 17 cm 2 but do not provide a refer- 
ence). 

That t coo \ ~ 4dv indicates our simulated column densities of hot 



4 For slowing down fast H in a sea of cold H + , there may also be a contri- 
bution to cr from "momentum transfer" in "elastic" (non-charge-exchange) 
collisions. This contribution increases c r over the viscosity cr oss section by 
only -30%; compare Figures 6 and 7 of[Schultz et al. 1 2008|). 



neutral hydrogen may be too large, but hopefully not by factors 
of more than a few. Keeping more careful track of the velocity 
distributions — and excitation states — of neutral hydrogen in the 
mixing layer would not only improve upon our calculations of 
Lyman-Qf absorption, but would also bear upon the recent detec- 
tion of B aime r Ha absorption in the hot Jupiters HD 209458b and 
HD 189733b (Ijensen et all2012n . 
(ii) Magnetic fields. Insofar as our results depend on Kelvin- 
Helmholtz mixing, our neglect of magnetic fields is worrisome 
because magnetic tension can suppress the Kelvin-Helmholtz in- 
stability dFranketal.il 199^) . For numerical simulations of magne- 
tized planeta r y wind s interacting with magnetized stellar winds, see 



tizea planetary winds in 
ICohen et all 120 11 al ibi). 



These magnetohydrodynamic simulations 
can track how planetary plasma is shaped by Lorentz forces, but 
as yet do not resolve how the planetary wind mixes and exchanges 
charge with the stellar wind. 

(iii) Dependence of Ly-a absorption A on the planetary wind density 
n p . In the same vein as item (ii), we found empirically that A oc n]J 2 , 
and argued that this result arose from the Kelvin-Helmholtz growth 
timescale. Ekenback et al. (2010) found a much weaker depen- 
dence: increasing n p by a factor of 100 only increases A in their 
models by a factor of ~2 at -100 km/s and even less at positive 
velocities — see their Figures 8 and 9. The true dependence of A on 
n p remains unclear. 

(iv) Rotational effects and gravity. There are a few order-unity geo- 
metrical corrections that our study is missing. Our standard stel- 
lar wind velocity of v* = 130 km/s is comparable to the planet's 
orbital velocity of v orb = 150 km/s, so that in reality the stellar 
wind strikes the planet at an angle of roughly 45 deg. The Cori- 
olis force will also deflect the planetary wind by an order-unity 
angle after a dynamical time of r/v orb ~ 5 x 10 4 s, by which time 
the wind will have travelled ~5R P from the planet. These geomet- 
rical e ffect s are potentially observ able — see, e.g., ISchneiter et al.l 
d2007h andlEhrenreich et all d2008l) for modeling of HD 209458b, 
and lRappaport et al] d2012h for a real-life example of a transit light 
curve that reflects the trailing comet-tail-like shape of the occulting 
cloud. However, these geometrical effects seem unlikely to change 
the basic order of magnitude of the absorption A ~ 10% that we 
have calculated. 

We have also neglected planetary gravity, stellar tidal gravity, and 
the centrifugal force, all of which can change the planetary wind 
velocity. But this omission seems minor, since we have drawn our 
input planetary wind velocities from calculations that do account 
for such forces (M09), at least along the substellar ray. According 
to Figure 9 of M09, the planetary wind accelerates from v p « 10 
km/s at a planetocentric distance d = 4R P , to v p « 30 km/s at 
d = \0R P . This range of velocities and corresponding distances 
overlap reasonably well with the range of velocities and distances 
characterizing our simulations. 

(v) Hydrodynamic approximations for the stellar and planetary 
winds. We have not formally justified our use of the hydrodynamic 
equations to describe the wind-wind interaction. The problem is 
that the collisional mean free path in the stellar wind is much longer 
than the lengthscales of the flow: /l C ouiomb,* = l/(^*crcouiomb) ~ 
10 13 (10 4 cm" 3 /^*)(10 -17 cm 2 /cr Co uiomb) cm, where cr Coulomb 
10 -17 (r*/10 6 K)~ 2 cm 2 is the cross section for protons scattering 
off protons. That the Solar wind is collisionless and does not nec- 
essarily admit a one-fluid treatment is well-known. 
Neve rtheless, it is perhaps just as well-known that Parker's dl958l 
Il963h use of the fluid equations to describe the collisionless 
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Solar wind is surprisingly accurate, capturing the leading-order 
features of the actual Solar wind. The role of Coulomb colli- 
sions in relaxing the velocity distribution functions of protons 
and electrons is fulfilled instead by plasma instabilities and wave- 
particle interactions — see , e.g., reviews of Solar wind physics by 
iMarsch et all d20Q3h and iMarschl d2Q06h . The gross properties of 
collisionless shocks can still be modeled with the hydrodynamic 
equations insofar as those properties depend only on the macro- 
scopic physics of mass, momentum, and energy conservation, and 
not on microphysics (e.g., Ishulll992h . 

Note that the planetary wind is fully collisional because 
of its higher density and lower temperature, and model- 
ing it as a single fluid appears justified: A Co u\omb,p ~ 
10 7 (10 6 cm -3 /n p )(T P /10 4 K) 2 cm, which is smaller than any other 
length scale in the problem. 



(vi) Non-Maxwellian behavior of the stellar proton velocity distribu- 
tion. Lyman-Qf absorption at the redshifted velocity of +100 km/s 
arises from charge-exchanged neutral hydrogen at the assumed stel- 
lar wind temperature of 10 6 K. We have assumed a Maxwellian 
distribution function for hydrogen in the stellar wind, and have 
ignored non-Maxwellian features that have been observed in the 
actual Solar wind, including high-energy tails and temperature 
anisotropics. Accounting for non-Maxwellian behavior may in- 
troduce order-unity corrections to our results for the absorption. 
For the more polar fast Solar wind, proton temperatures paral- 
lel to and perpendicular to the Solar wind magnetic field differ 
by factors of a few at heliocentric distances of 5-10 Solar radii 
(Mc Kenzie et al.lll997h . For the more equatorial slow Solar wind — 
which our simul ations are modeled a fter — temperature anisotropics 
are more muted (IMarsch et al.l l2003, page 391). 

(vii) Stellar radiation pressure. Stellar Lyman-a photons can ra- 
diatively accelerate neut ral hydrogen away from the star (e.g., 
IVidal-Madiar et alJliool M09). Both the planetary wind, and the 
charge-exchanged stellar wind in the mixing layer, are subject to a 
radiation pressure force that exceeds the force of stellar gravity by 
a factor p on the order of unity. 

Radiative repulsion of the charge-exchanged stellar wind in the 
mixing layer may not be observable, because once hot neutral hy- 
drogen is created in the mixing layer, it is advected off the projected 
limb of the star before radiation pressure can produce a significant 
velocity: Sv md ~ GM*/r 2 x p x t adv ~ 6p km/s, which does not 
exceed the hot neutral hydrogen's thermal velocity of ~100 km/s. 
What about radiative acceleration of the planetary wind? The travel 
time of the planetary wind from the planet to the mixing layer 
is ~10R p /v p ~ 10 5 s, long enough for neutral hydrogen to attain 
radiative blow-out velocities in excess of 100 km/s. However, the 
amount of hydrogen that suffers radiative blow-out is limited to the 
column that presents optical depth unity to Lyman-ar photons. This 
column is l/cr line _ ctr ~ 2 x 10 13 (r p /10 4 K) 1/2 cm -2 , and is much 
smaller than the typical column in the planetary wind, which is 
(1 - fp)n p R p ~ 10 16 cm -2 . Thus the bulk of the planetary wind 
is shielded from radiative blo w-out, and our neglect of radiation 
pressure appears safe. Note that lLecavelier des Etangs et aDd2012h 
find that radiation pressure cannot explain the largest blueshifted 
velocities observed for HD 189733b; like us, they favor charge ex- 
change. 
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